/*
 * F_gamma.cpp
 *
 *  Created on: 13/09/2012
 *      Author: rogsoares
 *
 *      For 2-D domains: get boundary edges contribution.
 */
#include "EBFV1_elliptic.h"
namespace PRS{
	int EBFV1_elliptic::gradient_F_gamma(Mat F, TMesh* pTMesh_data){
		static bool key = true;
		if (key){
			double Fij_column1[4], Fij_column2[4];
			int nelem, n, nedges;
			int dim = pTMesh_data->getMeshDim();
			int ndom = pTMesh_data->getNumDomains();
			for (int dom=0; dom<ndom; dom++){
				// OMEGA: all edges domains
				pTMesh_data->edge_list[dom].getSize(nedges,n);
				for (int i=0; i<nedges; i++){
					EdgeData* edata = pTMesh_data->edge_list[dom].getValue(i,0);				// get next domain's edge
					VertexData* vdata_0 = edata->vertex_0;										// edge's vertex 0
					VertexData* vdata_1 = edata->vertex_1;										// edge's vertex 1
					int ID_0 = vdata_0->ID_PETScOrdering;										// vertex 0 ID for matrix assembly
					int ID_1 = vdata_1->ID_PETScOrdering;										// vertex 1 ID for matrix assembly
					gradient_F_omega(pMatVec->F[dom],dim,ID_0,ID_1,edata);
				}
				// GAMMA: boundary domains
				pTMesh_data->bdry_elem_list[dom].getSize(nelem,n);
				if (dim==2){
					for (int i=0; i<nelem; i++){
						BdryElementData* edata = pTMesh_data->bdry_elem_list[dom].getValue(i,0);
						double *Dij = edata->Dij;
						double volumeI = edata->vertex_0->volume;
						double volumeJ = edata->vertex_1->volume;

						// fill edge matrix
						Fij_column1 [0] = 5.*Dij[0]/(6.*volumeI);
						Fij_column1 [1] = 5.*Dij[1]/(6.*volumeI);
						Fij_column1 [2] = Dij[0]/(6.*volumeJ);
						Fij_column1 [3] = Dij[1]/(6.*volumeJ);
						Fij_column2 [0] = Dij[0]/(6.*volumeI);
						Fij_column2 [1] = Dij[1]/(6.*volumeI);
						Fij_column2 [2] = 5.*Dij[0]/(6.*volumeJ);
						Fij_column2 [3] = 5.*Dij[1]/(6.*volumeJ);

						// where edge matrix must be assembled.
						int id0 = edata->vertex_0->ID_PETScOrdering;
						int id1 = edata->vertex_1->ID_PETScOrdering;
						int pos1 = 2*(id0-1);
						int pos2 = 2*(id1-1);
						int idxm[4] = {pos1,pos1+1,pos2,pos2+1};
						int idxn[2] = {id0-1,id1-1};
						ierr = MatSetValues(F,4,idxm,1,&idxn[0],Fij_column1,ADD_VALUES);
						ierr = MatSetValues(F,4,idxm,1,&idxn[1],Fij_column2,ADD_VALUES);
					}
				}
				else{
					throw Exception(__LINE__,__FILE__,"Fails!");
				}
				PETSc_assemblyMatrix(pMatVec->F[dom]);
			}
			key = false;
		}
	}
}
